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Charge carrier concentrations in acceptor-doped proton-conducting perovskites are to a large ex¬ 
tent determined by the hydration and oxidation of oxygen vacancies, which introduce protons and 
holes, respectively. First-principles modeling of these reactions involves calculation of formation 
energies of charged defects, which requires an accurate description of the band gap and the position 
of the band edges. Since density-functional theory (DFT) with local and semi-local exchange- 
correlation functionals (LDA and GGA) systematically fails to predict these quantities this can 
have serious implications on the modeling of defect reactions. In this study we investigate how the 
description of band gap and band edge positions affects the hydration and oxidation in acceptor- 
doped BaZrOs. First-principles calculations are performed in combination with thermodynamic 
modeling in order to obtain equilibrium charge carrier concentrations at different temperatures and 
partial pressures. Three different methods have been considered: DFT with both semi-local (PBE) 
and hybrid (PBEO) exchange-correlation functionals, and many-body perturbation theory within 
the GofFo-approximation. All three methods yield similar results for the hydration reaction, which 
are consistent with experimental findings. For the oxidation reaction, on the other hand, there is a 
qualitative difference. PBE predicts the reaction to be exothermic while the two others predict an 
endothermic behavior. Results from thermodynamic modeling are compared with available exper¬ 
imental data, such as enthalpies, concentrations and conductivities, and only the results obtained 
with PBEO and GofFo, with an endothermic oxidation behavior, give a satisfactory agreement with 
experiments. 

PACS numbers: 82.45.Un, 71.15.Mb, 71.20.Ps, 82.60.Cx 


I. INTRODUCTION 

Since the beginning of the 1980s, when Iwahara et alr^ 
discovered proton conduction in acceptor-doped SrCeOa, 
perovskite oxides {ABO 3 ) have been studied extensively 
with respect to their potential as proton conductors 
Such materials have many applications, including fuel 
cells, electrolyzers, hydrogen separation membranes and 
hydrogen sensors 4 Many acceptor-doped perovskites are 
also oxide ion, electron and hole conductors^^l and suit¬ 
able for applications such as electrodes and hydrogen sep¬ 
aration membranesi^ As some applications rely on the 
perovskite being a pure ionic or electrical conductor while 
others do not it becomes important to understand and 
control the conductivity mechanisms in order to predict 
and optimize the material performance. 

Proton incorporation into the perovskite structure is 
made possible through acceptor doping. By substituting 
S-site cations with dopant ions of lower valency posi¬ 
tively charged oxygen vacancies are formed due to charge 
compensation. By exposing the doped perovskite to wa¬ 
ter vapor the vacancies can be filled by water molecules 
which introduces protons into the structure. In Krdger- 
Vink notation this reaction is expressed as 

H20(g)+v^*+0g^20H^, (1) 

which describes how a water molecule, an oxygen va¬ 
cancy and an oxide ion form two hydroxide ions (pro¬ 
tons). Oxygen vacancies also enable the incorporation 
of holes, which can be introduced through oxidation of 


oxygen vacancies, 

1 02 (g)+v^*^2h* + Og. (2) 

Theoretical modeling based on density-functional the¬ 
ory (DFT) has become an important computational 
tool in materials science. The local density approxima¬ 
tion (LDA) and various semi-local generalized gradient 
approximations (GGAs) are routinely being used. In 
condensed matter research the Perdew-Burke-Ernzerho^ 
(PBE) type of GGA is currently the most common 
parametrizationi^ GGAs have been applied to study 
hydratioiii2r2£ as well as oxidatio n^"^d^~ — in different 
perovskite oxides. 

For the oxidation process the semi-local GGAs pre¬ 
dict the reaction to be exothermic d ^ dd "^^ 2 ^ the hole 
concentration is decreasing with increasing temperature. 
The hole conductivity is proportional to both the hole 
concentration and the hole mobility, and it is experi¬ 
mentally established that the hole conductivity increases 
with temperaturei^^^^ For the GGA result to be consis¬ 
tent with the experimental results it has therefore been 
suggested that the hole mobility increases more rapidly 
than the decrease in hole concentration]This, how¬ 
ever, is not in line with the common view in the research 
field and it has been stated that the electronic structure 
of acceptor-doped proton-conducting perovskites remain 
surprisingly poorly understoodi^ 

It is well-known that local and semi-local function¬ 
als underestimate the band gap of semiconductors and 
insulatorsa shortcoming that extends to the de¬ 
scription of the valence and conduction band edges. The 
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position of the top of the valence band is decisive for a 
correct description of the oxidation reaction in Eq. @ 
and hence one has to go beyond standard DFT with 
LDA/GGA in order to describe hole conduction prop¬ 
erly. 

In this paper we have performed theoretical modeling 
of hydration and oxidation of an acceptor-doped oxide. 
Three defects are of interest in this context, the doubly 
positively charged oxygen vacancy, the hole and the pro¬ 
ton, where the latter is often regarded as a hydroxide ion. 
Here we treat the hole as a band state. The method¬ 
ology is applied to acceptor-doped BaZrOs, one of the 
most promising proton conducting perovskites since it 
combines high bulk proton conductivity with chemical 

stability 

First-principles calculations are used to determine the 
electronic structure and defect formation energies. The 
starting point is DFT based on the PBE functional for 
the exchange-correlation energy. To remedy the band 
gap problem we then consider two different approaches. 

The first one is based on a many-body perturbation 
technique. We determine the quasiparticle correction to 
the PBE energy levels using the GoWo-approximation in¬ 
troduced by Hedin.— The second approach is based on 
hybrid functionals that admix a fraction of non-local ex¬ 
change to a semi-local approximation. We use the hybrid 
functional PBEO, which is obtained from PBE by replac¬ 
ing 25% of the PBE exchange energy by Hartree-Fock 
exchange. To make our study less empirical we stick to 
this original suggestion of 25% Hartree-Fock exchange^ 
and we do not make use of range-separation, as intro¬ 
duced in the corresponding HSE functionals^i Addition¬ 
ally, PBEO has been shown to give a good description of 
BaZrOs-^ 

Thermodynamic modeling based on our first-principles 
results is then performed to obtain charge carrier concen¬ 
trations in the acceptor-doped system at different tem¬ 
peratures and environmental conditions. We find a quali¬ 
tative difference for the oxidation reaction, being exother¬ 
mic with PBE and endothermic using the GqWq ap¬ 
proach and the hybrid functional. Indeed, only the latter 
behavior is found to be consistent with experimental data 
of charge carrier concentrations and hole conductivities. 

The paper is organized as follows. Section HIl describes 
the different aspects of the theoretical framework used in 
the paper while Section uni contains the computational 
details of the PBE, PBEO and Go Wo calculations. The 
results are presented and discussed in Section 113 and 
El and finally, a summary of the paper together with 
conclusions is given in Section IVII The Appendix gives a 
description of band structure alignment with respect to 
the vacuum level based on surface calculations. 


II. THEORY 

In this work we study the thermodynamics of defect 
configurations in the dilute limit. To this end, the for¬ 


mation free energies of individual point defects are calcu¬ 
lated (if necessary for different charged states) as a func¬ 
tion of atomic and electronic chemical potentials. The 
properties of the real system, most importantly defect 
concentrations, are then obtained by invoking the charge 
neutrality condition, which is employed to fix the elec¬ 
tronic chemical potential under different environmental 
conditions (atomic chemical potentials). An extensive 
review on the subject of first-principles modeling of de¬ 
fect formation in solids can be found in Ref. [3^ 

A. Defect formation energies 

The formation energy of a defect in charge state q is 
given by 

AEdef = + El,,, - ^ 

i 

+ 9(evBM + Me + An"^), (3) 

where E^\ and are the total energies of the defective 
and ideal systems, respectively. Am denotes the change 
in atomic species i upon defect formation and Mi is the 
corresponding chemical potential. Finally, fx, represents 
the electron chemical potential with respect to the va¬ 
lence band maximum, evBM- The terms E^,,, and Av‘‘ 
are corrections that compensate errors associated with 
charged defectsi^ The former term corrects errors due 
to image charge interactions, which are consequences of 
the periodic boundary conditions. The latter so-called 
potential alignment term corrects for the offset of elec¬ 
trostatic potentials of the charged defective and neutral 
ideal system. 

The band gap problem of DFT affects the formation 
energies and can be approximately corrected for by us¬ 
ing quasi-particle energy shifts from GqWo calculations. 
The method considered here, which is a perturbative ap¬ 
proach based on the DFT result, corresponds to applying 
a band gap correction to Eq. (|3]) and is described in more 
detail in Refs. In general, this approach requires 

knowledge of the shifts of both band edges as well as de¬ 
fect levels. Fully ionized defects, which is the nature of 
the defects in this paper, are only affected by the shift of 
the valence band edge. The band gap corrected forma¬ 
tion energy for such defects is given by 

= AESr + 9AevBM, (4) 

where AevBM = fvM - «VBM' 

For finite temperatures and pressures, Eq. © can be 
written as 

AGdef = G‘:‘f + El„, - G‘f - ^ Am9^ 

i 

+ g(evBM + Me + Az;'^), (5) 

where Gjj°f and Gf^* are the Gibbs free energies of the 
defective and ideal systems respectively, and gt is the 
chemical potential of the elemental reference phase i at 
finite temperatures and pressures. 
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B. Chemical potentials of the gas phase 

The considered defects are oxygen vacancies and pro¬ 
tons and chemical potentials for O and H are therefore 
needed. The environments of the oxidation and hydra¬ 
tion reaction are oxygen gas (O 2 ) and water vapor (H 2 O), 
and the chemical potentials of O and H are thereby ex¬ 
pressed as 


where is the electronic contribution, and the tem¬ 
perature dependent terms and S''^'^{T) represent 

vibrational contributions. The latter two are calculated 
within the harmonic approximation using an Einstein 
modelJ^di Here we assume that the formation of a defect 
does not affect the vibrational frequencies of neighboring 
atoms. The change in and S^'^{T) due to the 

addition of one atom of species i is given by 


go = ^902 (6) 

5H = 25 H 2 O - ^502- (7) 

By assuming an ideal gas behavior the chemical potential 

of O 2 at temperature T and partial pressure P 02 (and 
equivalently for H 2 O at PH 20 ) can be written as 

g02(T,P02) = 902 + + ^ 02 (^) - Ts°o,{T) 

+ kT\n^, (8) 

P 02 

where SqJ' is the zero-point energy of the O 2 molecule 
and (T) and Sq^ (T) represent the temperature depen¬ 
dencies of enthalpy and entropy of the gas phase at the 
reference pressure Pq^ . The enthalpies and entropies of 
O 2 and H 2 O are extracted from thermodynamic tablesi^ 
Within the harmonic approximation the zero-point ener¬ 
gies are given by where Wfc are the molec¬ 

ular vibrational frequencies. Experimentally determined 
frequencies^!^ yield Sq^' = 0.10 eV and = 0.56 eV. 

Total energies from DFT are used for Common 
practice is to use the molecular total energies 

902 = 7^02 (9) 

/^H20 =-^H2*0' (10) 

This is problematic since PBE is known to overbind the 
O 2 molecule with 0.9 eV. To overcome this problem total 
energies of atoms are used instead and combined with ex¬ 
perimental values for the cohesive energies according 
to 

P02=2El°^+eh°2 ( 11 ) 

MH 20 = 2 i^H‘ + l?o‘+^H 20 - ( 12 ) 

With experimental data from Ref. we obtain = 

—5.21eV and £h 2 ^ = —10.07eV, where the zero-point 
energies (see above) have been removed. 


C. Free energy of the solid phase 


AUf'^iT) = ^ 

k^l 

ASf\T) = kJ 2 






1 2 gfi(^i,klkT _ ^ 


(14) 


/C=l ^ 


huJi,k/kT / 

-huji k/kT \ 

(,huJi,klkT _l 

1 e JJ 


(15) 


where uji^k are the vibrational frequencies. For the oxy¬ 
gen atom we use the frequencies 557cm“^, 250 cm“^ and 
250 cm“^, and for the proton we use 3502 cm“^, 900 cm“^ 
and 601 cm“^, which have been extracted from Ref. M 
andflTl 


D. Defect concentration 


Defect concentrations are considered to be within the 
dilute limit and are therefore given by 


Cdef — 


Ndef 


„-AGdef/feT 

^ 5 


(16) 


where iVdef is the number of defect sites in the primitive 
cell with volume Vq. In this case Vc = Og with ag being 
the lattice constant. There are three oxygen sites in the 
primitive cell and therefore three available sites for the 
oxygen vacancy, i.e., JVy = 3. Proton sites are associated 
with oxygen ions, with four possible configurations per 
oxygen site,— which yields A^h = 12 proton sites in each 
primitive cell. 

In order for the dilute-limit approximation to be valid 
the occupancy has to be much smaller than the num¬ 
ber of available sites (cdefK IVdef)- In this paper we 
use a dopant concentration of 10%, which yields a max¬ 
imum proton occupancy of 0.1 per primitive cell. This 
corresponds to 1 in 120 proton sites being occupied. The 
same dopant concentration yields an upper limit of 0.05 
oxygen vacancies per primitive cell, which corresponds 
to 1 in 60 oxygen sites being vacant. The dilute-limit 
approximation is thus justified. 


The considered expression for the free energy of the 
solid phase depends only on temperature since the PE- 
term is very small within this context and can be ne¬ 
glected. This implies that the Gibbs and Helmholtz free 
energies are practically identical and one can write 

G(T) « F(T) = -f t7''''"(r) - r5'"'"(T), (13) 


E. Electron chemical potential 

The electron chemical potential pe is obtained by solv¬ 
ing the charge neutrality condition 

^ kl^i^Po) -f TI\i(^Pq^ — 0, 

def 


( 17 ) 
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where and are the electron and hole concentration, 
respectively, and the sum is over all defects in the mate¬ 
rial including the acceptor dopants. Equation (El) can 
be solved by iteration4i In the present work electrons 
and holes are treated as band states and the correspond¬ 
ing concentrations are obtained from the density of states 
(DOS) 5 (e) according to 


rie = / g(e)/(e,/ie)de (18) 

eCBM 

/ eVBM 

5 (e) [1 -/(OMe)] rfe, (19) 

-00 


where evBM and ccbm denote the positions of 
the valence band maximum (VBM) and conduction 
band minimum (CBM), respectively, and /(e,/Xe) = 
{exp [(e — evBM — f^e)/kT] + 1}“^ is the Fermi-Dirac dis¬ 
tribution function. The DOS is determined from first- 
principles calculations. 
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Figure 1. Convergence of quasi-particle energies from GqWo 
calculations based on PBE wave functions with respect to the 
number of bands included in the calculation. 
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III. COMPUTATIONAL DETAILS 

First-principles calculations within the density- 
functional theory (DFT) formalism were carried out 
using the Vienna ab-initio simulation package;^ which 
uses plane wave basis sets and periodic boundary condi¬ 
tions. The projector augmented wave method (PAW)^ 
was employed to describe ion-electron interactions. 
Two different functionals were used to model exchange 
and correlation in their non-spin polarized versions: 
the generalized gradient approximation functional 
PBE^ and the hybrid functional PBEOi^ The plane 
wave cutoff energy was set to 400 eV and a 6 x 6 x 6 
Monkhorst-Pack grid was used for /c-point sampling of 
the BaZrOa primitive cell and then reduced accordingly 
with increasing supercell size. Super cells comprising up 
to 6 X 6 X 6 unit cells were used for defect calculations 
based on the PBE functional. PBEO calculations were 
conducted for 3x3x3 supercells only. Ionic relaxation 
was carried for all structures until the residual forces 
were below 0 . 02 eVA“^. 

All calculations were performed with the cubic per- 
ovskite structure, which belongs to space group Pm3m. 
The optimized PBE lattice constant of 4.235 A is some¬ 
what larger than experimental values 4.191-4.197 Aj^i^ 
but in agreement with previous theoretical studies of 
BaZrOa based on GGA functionalsj ^^d^d'^ The PBEO cal¬ 
culations were carried out at the PBE lattice constant for 
a more direct comparison. 

Many-body calculations were carried out within the 
formalism of the quasi-particle method GW ^ More 
specifically, the Go Wo approach was used. Galculations 
were based on PBE wave functions and employed PAW 
data sets optimized for GW calculations.— The general 
plane wave cutoff energy was 434 eV while a cutoff of 
290 eV was employed in the response function calcula¬ 
tions. The Brillouin zone was sampled using a P-centered 


5x5x5 fc-point mesh and all calculations were carried 
out at the PBE lattice constant. 

While the band gap converges relatively quickly with 
the number of empty states included in the calculations, 
individual quasi-particle energies typically converge more 
slowly. As shown in Fig. [1] VBM and CBM are, how¬ 
ever, observed to depend linearly on the inverse number 
of bands, whence converged values were obtained by ex¬ 
trapolation. This a ppr oach is similar to the hyperbolic 
fit employed in Ref. 


IV. RESULTS 

A. Electronic structure 

The band structure of BaZrOa from PBE and PBEO 
calculations is presented in Fig. [5] The band gap is in¬ 
direct with the VBM at R and the CBM at P. The 
size of the gap, which is determined from single-particle 
eigenvalues;^ is 3.13 eV and 5.35 eV with PBE and 
PBEO, respectively. The direct band gap, with the VBM 
and CBM at P, is only slightly larger: 3.38 eV with PBE 
and 5.57eV with PBEO. The shape of the band struc¬ 
tures is very similar in both cases, which indicates that 
main difference between PBE and PBEO lies in the size 
of band gap and the position of band edges. This is illus¬ 
trated in Fig. [31 which shows the total (DOS) and partial 
density of states (PDOS). The valence band consists of 
oxygen p-states while the conduction band of zirconium 
d-states. 

To compare the position of band edges the PBE and 
PBEO band structures need to be properly aligned. Such 
an alignment can be done with respect to common ref¬ 
erence potential, e.g., the average local electrostatic po¬ 
tential or the vacuum level— In this study we use the 
same pseudopotentials and lattice constant for both PBE 
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PBE PBEO 




Figure 2. Comparison of PBE and PBEO band structures 
for BaZrOa. Blue and red lines represent empty and occupied 
bands, respectively. The grey areas indicate the extent of the 
indirect band gap (R-F). The energy scale is chosen to be 
zero at the PBE VBM. 
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Figure 3. Total and partial density of states for BaZrOs 
obtained with PBE and PBEO. The red and blue lines cor¬ 
respond to electronic p- and d-states and the dashed vertical 
line marks the VBM. 


and PBEO and the ionic contribution to the electrostatic 
potential is therefore the same. The electron density is 
found to be very similar with both methods, which yield 
similar contributions to the potential as well. As a con¬ 
sequence the average local electrostatic potential is ap¬ 
proximately the same for both methods and the two band 
structures should be properly aligned4S We have also 
performed alignment with respect to the vacuum level 
using surface calculations, which verifies this alignment 
(see Appendix). 

Band gaps and band edge positions are summarized in 
Table H] and visualized in Fig. 01 where the band edge po¬ 
sitions are given with respect to the PBE VBM. Both 
PBE-bGoVFo and PBEO open up the band gap, from 
3.13 eV to 4.73 eV and 5.35 eV, respectively, and yield 
VBM/CBM shifts that are qualitatively similar. The 
rather good agreement between PBEO and PBE-l-GolTo 
calculations for the VBM offset is not trivial as it has 


Table I. Comparison of theoretical and experimental band 
gaps Agap, as well as VBM and CBM shifts Ae obtained from 
PBEO and PBE-l-GoVFo calculations with respect to PBE cal¬ 
culations. All values are given in units of eV. The theoretical 
data are also visualized in Fig. 01 


Method 

AevBM 

AecBM 

Agap 

PBE 



3.13 

PBE-bGoWo 

-1.10 

0.50 

4.73 

PBEO 

-1.42 

0.80 

5.35 

Experiment 
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Figure 4. Size and relative position of band gap for PBE, 
PBEO and PBE-l-GoWo calculations, where zero is set at the 
PBE VBM. The dashed lines indicate that PBE-fGolFo is a 
perturbative approach based on PBE. Also see Table U 


been shown that PBEO band edge positions can differ 
quite substantially from GqWo calculations, especially 
for wide band gap materials In general one should ex¬ 
pect PBE-|-GotTo calculations to be more reliable for this 
purpose as they represent a more rigorous theoretical ap¬ 
proach. 

There are several experimental values for the band gap 
of BaZrOa in the literature. Robertson^ reports a value 
of 5.3 eV, which in close agreement with the PBEO re¬ 
sult. More recent studies by Cavalcante et al^ and 
Yuan et al^ report band gaps in the range 4.8-4.9eV, 
which agree better with PBE-l-GoWo. The fact that the 
PBE-|-GoWo still slightly underestimates the experimen¬ 
tal band gap is consistent with calculations on other wide 
band gap materialsi ^^’^°’^^ 


B. Defect formation energies 

Formation energies have been calculated for the oxy¬ 
gen vacancy AE^ and the proton AEr. The considered 
charge state of the vacancy is 4-2, which is the relevant 
state for the oxidation and hydration reactions. 

The terms and Av’^ in the expression for the for¬ 
mation energy (see Eq. ®) are corrections to errors in¬ 
troduced by charged defects and periodic boundary con- 
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Figure 5. Convergence of formation energies for the +2 
charged oxygen vacancy (v) and the proton (H) with respect 
to supercell size (number of atoms). Solid lines are fits of the 
data to Eq. The calculations are performed with the 

PBE functional. 


ditions. Several correction schemes have been proposed 
over the years to reduce these errors (see Refs. [H andl^ 
for examples). Here we employ the finite-size scaling ap¬ 
proach, in which the formation energy is calculated for 
several different supercell sizes and the corrected value 
Eoo is obtained by fitting the data points to a polyno¬ 
mial of the form 

E{N) = aN-^+bN-^^^+ Eoo, ( 20 ) 

where N is the number of atoms in the supercell. In this 
fashion not only the leading terms of the multipole expan¬ 
sion of the electrostatic image interaction^ are accounted 
for but also elastic image interactions)^ This approach is 
suitable in this case since it is computationally feasible to 
obtain a sufficiently large number of data points for a rea¬ 
sonable fit. Additionally, since the screening in BaZrOa 
is quite large (the static dielectric constant Sr has been 
experimentally measured to fall in the range 40"160^), 
electrostatic image charge interactions, which are propor¬ 
tional to e“^, can be expected to be small. There is thus 
no benefit in using more advanced schemes. 

Supercells with up to 6 x 6 x 6 unit cells are considered 
for the extrapolation, which corresponds to 1080 atoms 
in the non-defective configuration. The results for the 
PBE functional are shown in Fig. [5j The extrapolated 
formation energy for the oxygen vacancy is 1.31 eV while 
a value of 0.25 eV is obtained for the proton. The figure 
shows that the formation energy of both defects is quite 
close to the extrapolated value already for 3 x 3 x 3 su¬ 
percells (135 ± 1 atoms), which is related to the strong 
electrostatic screening. Since the PBE formation ener¬ 
gies of the 3x3x3 supercell are already very close to 
the extrapolated value, this supercell size was employed 
for PBEO calculations, which are computationally much 
more demanding. 

Defect formation energies obtained from PBE and 


Table II. Comparison of formation {AEv and AEu) and reac¬ 
tion {AEox and AEhydr) energies, where the former are given 
for the electron chemical potential being located at the VBM 
(/ie = 0). For PBE and PBE-|-x[GolFol the formation ener¬ 
gies are extrapolated values, see Eq. (1201) and Fig. [S] PBEO 
values correspond to 3 x 3 x 3 supercells. All energies are 
given in units of eV. 


Method 

AEx 

AAh 

^-E'ox 

^-^hydr 

PBE 

1.31 

0.25 

-1.31 

-0.82 

PBE-bx[GolFol 

-0.88 

-0.85 

0.88 

-0.82 

PBEO 

-1.53 

-1.22 

1.53 

-0.90 


PBEO calculations are summarized in Table HU All val¬ 
ues are determined at the VBM corresponding to = 0. 
The differences between the PBEO and PBE values are 
—2.84eV for the vacancy and —1.47eV for the proton. 
These differences are very close to 2AevBM and Acvbm 
(see Table m, which indicates that the difference between 
PBE and PBEO is mostly due to the shift of the VBM. 
This observation in turn validates the PBE-|-x[GoWo] ap¬ 
proach. 


C. Reaction enthalpies and entropies 

The energy of the oxidation reaction in Eq. m is de¬ 
termined according to 

AAox = 2^e - AEv(Aie), (21) 

which is independent of fie- Calculated values for AAqx 
are listed in Table im With PBE the oxidation energy is 
— 1.31 eV, which implies an exothermic reaction favoring 
the formation of holes. With PBE-|-x[GoWo] and PBEO 
the oxidation energy is 0.88 eV and 1.53 eV, respectively, 
which corresponds to an endothermic reaction favoring 
oxygen vacancy formation. 

The energy of the hydration reaction in Eq. o is given 
by 

AEhydr = 2AEH(/re) — AE^{fj,e), (22) 

which, like the oxidation energy, is independent of fie- 
All three methods predict the reaction to be exother¬ 
mic with a similar magnitude for AEhydr, see Table HIl 
The reaction is slightly more energetically favorable with 
PBEO compared to PBE, while PBE and PBE-|-x[GolTo] 
yield identical values by construction. This close agree¬ 
ment between the different methods can be traced to the 
fact that the hydration energy does not depend on the 
position of the VBM. 

The standard enthalpy for both reactions can be de¬ 
termined from AEok and AAhydr by including the zero- 
point energies and the temperature dependence of both 
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Figure 6. Temperature dependence of the standard enthalpy 
and entropy of the hydration and oxidation reactions, see 
Equations o and (H). The electronic contributions to the 
enthalpy (A_Eox and AShydr) have been subtracted, thus the 
values at zero temperature correspond to zero point energies. 


the solid and the gas phase. The enthalpies are given by 
AH:,{T) = + AC/J'^T) - iegf - 

(23) 

= Ai^hydr + 2AC/"^T) + AUS'^iT) 

~ ~ ^H2o(^)- (24) 

Similarly, the entropies are given by 

AS:^{T) = AS^S^T) - (T) (25) 

AS^^^,{T) = 2AS^\T) + AS^S\T) - s°^^oiT). (26) 

In Fig.[n]we show the standard enthalpy and entropy as 
a function of temperature for both reactions. AEox and 
Aifhydr have been subtracted from the enthalpy thus the 
values at zero temperature correspond to the net zero- 
point energy of the reactions. These values are much 
less than the zero-point energy of the respective phases, 
which indicate that there is a large cancellation effect. 
Thus, if zero-point motion effects are included it is of 
importance to consider contributions from both the gas 
and solid phases. 


D. Oxidation 

Based on the computed formation energies the equilib¬ 
rium defect concentrations can be determined for differ¬ 
ent temperatures and pressures using the self-consistent 
scheme described in SectionlTTl With these concentrations 


the oxidation reaction can be studied by calculating the 
corresponding equilibrium constant 


ifox(T) 


\ -1/2 2 
PO 2 \ 

PO2 ) Cv 


(27) 


where Cv and cq denote oxygen vacancy and oxygen ion 
concentrations, respectively. 

In Fig. [ 7 ] we show the equilibrium constant as function 
of temperature together with the hole concentration of 
a 10% acceptor-doped system at the reference pressure 
{PO 2 = Ibar). As can be expected from the oxidation 
enthalpies, the results differ quite significantly between 
PBE and the other two methods. With PBE the hole 
concentration increases with decreasing temperature and 
is completely compensating the dopant charge at lower 
temperatures. With PBEO and PBE-|-y[G'oWo] the con¬ 
centration displays the inverse temperature dependence 
and is several orders of magnitude smaller. These fea¬ 
tures are reflected in the equilibrium constant, where the 
positive slope of the PBE curve indicates an exothermic 
process while the negative slope obtained using the other 
two methods corresponds to an endothermic reaction. 

In general, the slope of the In iT(r)-curve is considered 
to correspond to the enthalpy of the reaction. We define 
an effective oxidation enthalpy according to 


ah:^^%t) = -k 


dlnifox(T) 

d{l/T) 


(28) 


Fitting the data in Fig.[7]to Eq. (1^ yields AE[°^^{T = 
1000 K) values of -0.66 eV, 1.30 eV and 1.92 eV for PBE, 
PBE-|-x[GoWo] and PBEO, respectively. These values 
can be compared with —1.26eV, 0.93 eV and 1.58 eV for 
AH°^{T) at T= 1000K. 

The electron chemical potential, which is also depicted 
in Fig. [71 is negative with PBE below 1000 K and re¬ 
mains close to the valence band edge for larger tempera¬ 
tures. For PBEO and PBE-|-x[GolFo] on the other hand 
the electron chemical potential is located well within the 
band gap over the entire temperature range. In the latter 
case the Boltzmann approximation can be used to find 
a more simplified expression for A'ox(T') and nu. The 
equilibrium constant can then be written as^ 

K,4T) = [nyB{T)f e-^^ox(T)/fcT^AS°.(T)/ic^ (29) 

where nvB(T) = 2(mj(fcT/27r?i^)^/^ and is the effec¬ 
tive mass for the hole. From this expression it follows^ 
that 


Ai7°i"®(T) = AH°^{T) + 3kT. (30) 

The contribution 3feT stems from the holes and is equal 
to 0.26 eV at 1000 K. This explains the difference be¬ 
tween the slopes of the PBE+x[GoWo] and PBEO curves 
in Fig. [7] and the corresponding oxidation enthalpies 
AH°^(T). While for the PBE there is also a positive 
contribution to AH°^{T) it is more difficult to obtain an 
explicit expressioni^ 
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Figure 7. The upper panel shows the equilibrium constant 
of the oxidation reaction in Eq. m while the middle and 
lower panels show the corresponding hole concentration and 
self-consistently obtained electron chemical potential. The 
concentrations are obtained with po^ = 1 bar and a dopant 
concentration of 10%, where the latter is depicted as a black 
dashed line in the middle panel. 


We have also studied the dry system for a wide range 
temperatures and oxygen partial pressures. In Fig. |8]we 
show the hole concentration for different temperatures 
and oxygen partial pressures at a dopant concentration 
of 10%. The holes completely compensate the accep¬ 
tor dopants at high partial pressures if PBE energies are 
used, and the hole concentration is still quite substan¬ 
tial when the pressure decreases. With PBE-|-x[G'oWo] 
and PBEO we obtain a different picture. Here does the 
hole concentration become large only at high tempera¬ 
tures and very high partial pressures, and consequently 
the acceptor dopants are compensated by oxygen vacan¬ 
cies over most of the considered range. 


E. Hydration 

In the same manner as for the oxidation reaction, the 
hydration reaction can be studied through the corre¬ 


sponding equilibrium constant 

PHsO 


KhydriT) = 


-1 




CyCO 


( 31 ) 


where ch is the proton concentration. In this case the 
equilibrium constant can be written as^ 


The difference in the number of sites available for protons 
(IVh) and oxygen vacancies (W) introduces an additional 
conhgurational contribution to the entrop} d^i^° and we 
can dehne an effective hydration entropy according to 

= + . (33) 


In the present case we have iVn = 4:Nv and the additional 
term is equal to 0.24meV/K. 

At 900 K we obtain hydration enthalpies of —0.68eV 
with PBE and PBE-|-x[GoWo], and —0.76eV with PBEO. 
The corresponding effective hydration entropy at the 
same temperature is —1.38meV/K. 


F. Experimental conditions 

The environmental conditions in experimental studies 
are often such that both hydration and oxidation take 
place simultaneously. This is the case for a hydrated 
material under oxidizing conditions and during such cir¬ 
cumstances it is not possible to consider the two reactions 
independently. 

We have employed the scheme described in Section HIl 
to model these experimental conditions. Concentration 
profiles for a 10% doped material under wet conditions 
with PH 2 O = 0-02 bar and P 02 = 10“® bar are shown in 
Fig. [9l The material is hydrated at lower temperatures 
according to all three methods but only completely pro- 
tonated for PBEO and PBE-l-GoHo. With PBE the hy¬ 
dration occurs in competition with hole formation lead¬ 
ing to a situation with roughly 50% protons and 50% 
holes. Similar to dry conditions, the hole concentration 
increases with increasing temperature for both PBEO and 
PBE-|-x[GoWo] while the behavior is the opposite for 
PBE. 

In this study only isolated defects are considered, 
which is reasonable for low dopant concentrations. How¬ 
ever, at higher concentrations defect ordering and associ¬ 
ation effects cannot be neglected. Real systems are often 
subject to high dopant concentrations of approximately 
20% and above. While in such situations defect-defect 
interactions should be included we do not consider this 
complication in the present work. The scheme employed 
here (Section |n| can, however, be extended in straight¬ 
forward fashion to account for additional defect species, 
including defect pairs, as a first order approximation to 
defect-defect interactions. 
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Figure 8. Hole concentration calculated based on PBE, PBEO and PBE+x[GolFo] data under dry conditions at different 
temperatures and oxygen partial pressures. The dopant concentration is 10%, which corresponds to 1.3 x 10^^ cm“®. 
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Figure 9. Concentration profiles calculated based on PBE, PBEO and PBE+x[GoHo] data under hydrated conditions at 
PH 2 O = 0.02 bar and po^ = 10“® bar. The dopant concentration is 10%, which corresponds to 1.3 x 10^^ cm“^. 


V. DISCUSSION 
A. Hydration 

It was established in the previous section that the three 
methods considered in this work all predict very similar 
results for the hydration reaction. This shows that a 
change in description of the electronic structure has a 
small effect on the hydration enthalpy. 

The hydration of acceptor-doped BaZrOa has been 
studied extensively by several experimental groups and a 
compilation of their results is provided in Table lllll There 
is good agreement between the different experimental re¬ 
sults obtained at higher temperatures with some slight 
differences due to doping, where the dopant species ap¬ 
pears to have a more prominent impact on the results 
compared to the dopant concentration. 


At T = 900 K the calculated hydration enthalpy is 
= -0.68 eV with PBE and PBE-bx[GotTo], and 
= —0.76eV with PBEO. Since these values are 
computed for an effectively acceptor-doped BaZrOa sys¬ 
tem there is no specific entry in Table imi to compare 
with, although the values do agree quite well in gen¬ 
eral. For the same temperature the calculated effective 
hydration entropy is = — 1.38meV/K. The mag¬ 

nitude of this value is somewhat larger than the exper¬ 
imental entropies listed in Table uni Recent investiga¬ 
tions have shown that a more accurate treatment of the 
lattice vibrations gives a considerably better agreement 
with experimentsi^ 

There is one entry in Table uni which differs from the 
others, namely the 20% yttrium-doped system studied 
at low temperatures by Yamazaki et al.^ The absolute 
value of the enthalpy is much smaller in this case, which 
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Table III. Experimental values of hydration enthalpies and entropies for various acceptor-doped BaZrOs systems. 


System 

T (K) 

PH 2 O (bar) 

AHhAdr (eV) 

ASOxdr (meVK-i) 

Reference 

BaZro. 98 Yo. 0203-5 

773-1173 

0.023 

-0.84 

-0.98 

Kreuer et al^ 

BaZro.gsYo.osOa-a 

773-1173 

0.023 

-0.82 

-0.97 

Kreuer et al^ 

BaZro.gYo.iOs-i 

773-1173 

0.023 

-0.82 

-0.92 

Kreuer et aZ.— 


773-1073 

0.005-0.04 

-0.77 ±0.03 

-0.90 ±0.10 

Schober & Bohn^ 


573-1173 

O.l-l.O 

-0.84 ±0.04 


Kjqlseth et al^ 


673-873 

0.02 

-0.86 

-0.95 

Ricote et al^ 

BaZro.ssYo.isOs-i 

773-1173 

0.023 

-0.86 

-0.95 

Kreuer et al^ 

BaZro.sYo.gOs-i 

773-1173 

0.023 

-0.97 

-1.07 

Kreuer et aZ.— 


323-773 

0.023 

-0.23 ±0.01 

-0.40 ±0.01 

Yamazaki et al^ 


773-1173 

0.023 

-0.73 

-1.04 

Yamazaki et al^ 

BaZro.gSco.iOs-a 

773-1173 

0.023 

-1.24 

-1.29 

Kreuer et al^ 

BaZro.gGdo.iOa-^ 

773-1173 

0.023 

-0.69 

-0.89 

Kreuer et aZ.— 

BaZro.gluo.iOg-i 

773-1173 

0.023 

-0.69 

-0.93 

Kreuer et aZ.— 


corresponds to a less exothermic reaction. The authors 
argue that the difference with respect to other results is 
that the hole concentration can be neglected at low but 
not at high temperatures. This explanation is not con¬ 
sistent with either the PBE or the PBE0/PBE-|-x[GoWo] 
results in Fig. [HI Kjplseth et al^ on the other hand argue 
that the less exothermic behavior is due to association 
and ordering between defects and dopants, under the as¬ 
sumption that oxygen vacancies are more associated and 
ordered compared to protons. 


B. Oxidation 

The modeling of the oxidation reaction yields very dif¬ 
ferent results depending of the method that is consid¬ 
ered. The values of the oxidation enthalpy in Table HIl 
show that the standard DFT approach based on the PBE 
exchange-correlation functional predicts the reaction to 
be exothermic while PBE-|-x[GolEo] and PBEO predict 
an endothermic behavior. 

Figures [5] and [5] show that the exothermic nature of 
the PBE results yields large hole concentrations. The 
results for PBE in the latter figure indicate that 50% of 
the oxygen vacancies are oxidized even under hydrated 
conditions. This is inconsistent with experiments, where 
almost completely hydrated samples are obtainedi^*^ 

Unlike for the hydration reaction, there are to our 
knowledge no reported experimental values of the oxida¬ 
tion enthalpy for BaZrOa systems in the literature. There 
are however experimental values for other perovskite ox¬ 
ides, namely BaCeOa^ BaTiOg^ and SrTiOai^ The 
oxidation enthalpies for these systems (see Table HV)) are 
all positive, which corresponds to the reaction being en¬ 
dothermic. 

To compare these experimental values with theoreti¬ 
cal predictions AEqx was calculated for these perovskites 
as well. Calculations were performed with both PBE 
and PBEO using the same computational setup as for 
BaZrOa. Band gaps and band edge shifts were deter¬ 


mined as well, where the latter were obtained under the 
assumption that the PBE and PBEO band structures are 
aligned. Although the cubic perovskite structure is not 
the ground state for these materials it was chosen for 
simplicity. 

The results of these calculations are shown in Table HVl 
These three perovskites behave qualitatively similar to 
BaZrOa with negative and positive oxidation energies 
with PBE and PBEO, respectively. The latter are in bet¬ 
ter agreement with the experimental data. The band 
gaps are also improved for all systems and the VBM and 
CBM are shifted downwards and upwards respectively 
for all materials, similar to BaZrOa. The fact that the 
overall improvement of PBEO over PBE is a general fea¬ 
ture for these three systems in combination with their 
similarities to BaZrOa strongly suggests an endothermic 
oxidation reaction in BaZrOa. Thus, going beyond stan¬ 
dard DFT is a necessary procedure when studying the 
oxidation reaction in these materials. 

Throughout this article we have considered the hole to 
be a delocalized band state. If the hole instead would 
be a localized polaronic state (small polaron) then the 
oxidation enthalpy would be reduced by the formation 
energy of the polaron. Recent theoretical studie a^^i^^ 
based on the HSE functional and LDA-I-C/ show indeed 
that polaron formation is favorable in several perovskites 
(SrTiOa, BaTiOa and CaTiOa). However, the polaron 
formation energies are only about 0.1-0.2eV and thus 
quite small. While polaron formation would reduce AAox 
by 0.2-0.4eV it would not change the main conclusions 
of the paper. 

C. Conductivity 

Conductivity is a quantity that can be experimentally 
measured much more easily than defect concentrations. 
The conductivity of a charge carrier i can be decomposed 
into 




( 34 ) 
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Table IV. Lattice constants ao, band gaps Sgap, band edge shifts and oxidation enthalpies AEox /for several perovskite 
oxides. All calculations have been performed with the cubic perovskite structure. For the band edge shifts it is assumed that 
the PBE and PBEO band structures are aligned. Energies and lattice constants are given in units of eV and A, respectively. 


System 

ao 


Egap 


Band-edge shifts 

A£/ox 


PBE 

Exp. 

PBE 

PBEO 

Exp. 

VBM 

CBM 

PBE 

PBEO 

Exp. 

BaCeOs 

4.476 

4.4452^ 

2.25 

4.95 

4.41-21 

-1.45 

1.25 

-1.35 

1.85 

TJ^ 

BaTiOs 

4.031 

3.991-2. 

1.71 

3.82 

3.2& 

-1.41 

0.70 

-1.07 

1.96 

0.922L 

SrTiOs 

3.939 

3.9002^ 

1.81 

3.98 

3.2522 

-1.41 

0.77 

-1.00 

1.87 

1.40^ 

BaZrOs 

4.235 

4.191-^ 

3.13 

5.35 

4.8-5. 32.-.^^ 

-1.42 

0.80 

-1.31 

1.53 



where qi is the carrier charge, Bi the mobility and rii is 
the carrier concentration. 

Total and partial conductivities of yttrium-doped 
BaZrOs have been determined experimentally by several 
research groups.— With a fit to the Arrhenius like 
expression 

Tuh = (35) 

the reported hole conductivities Uh yield activation ener¬ 
gies Ea in the range 0.62 eV to 1.05 eV—^— To compare 
the experimental hole conductivities with our results for 
the hole concentrations the mobility of the holes is re¬ 
quired. While the mobility and hence the diffusion coef¬ 
ficient have been experimentally determined for both pro¬ 
tons and oxygen vacancies in yttrium-doped BaZrOa, the 
hole mobility is unknown. There are, however, mobilities 
reported in the literature for other perovskites including 
BaTiOg^ and SrTiOa— In Fig. [TU] these hole mobili¬ 
ties are depicted together with the proton and oxygen 
ion mobility in BaZro.gYo.iOa-^ based on experimental 
data from Kreuer et oL— Unlike the proton and oxygen 
ion mobilities, which clearly show temperature activated 
behavior, the hole mobilities have a temperature depen¬ 
dence close to T~^ corresponding to scattering limited 
band conduction mechanism. 

By assuming that ^ T~^ it follows from Eq. (1341) 
and Eq. (1351) that rih ^-Ea/kT^ consider 

PBE-|-x[GoWo] and PBEO, where rih Cv and the Boltz¬ 
mann approximation is valid, we get^ 


r(K) 


1500 1000 700 500 400 



iooo/r(K-i) 

Figure 10. Experimental mobility of charge carriers. The 
proton and oxide ion mobility is for BaZro.gYo.iOs-a and is 
based on data from Ref. [^ . The hole mobilities are based 
on the expressions given in Ref. (BaTiOs) and Ref. 
(SrTiOa). 


activated process involving small polarons with a mobil¬ 
ity given by B^ ^ T~^ exp{—E^ig/kT). In the present 
case the activation energy for hole migration E^ig has to 
be at least 1 eV, which is unlikely. 


VI. SUMMARY AND CONCLUSIONS 


At T = 1000 K the calculated oxidation enthalpies yield 
Ea = 0.65 eV and Ea = 0.96 eV for PBE-hx[GoWo] and 
PBEO, respectively, which are within the range of the 
experimental results— “— 

On the other hand, if we consider PBE the oxidation 
reaction is exothermic and Ea is negative (c.f. Fig. 0. 
This can not be made consistent with the measured con¬ 
ductivity under the assumption of a weakly temperature 
dependent mobility, B^ ~ T~^. For the PBE result to 
become consistent one has to assume a str ong ly temper¬ 
ature dependent mobility. In Refs. [13 andliol it was sug¬ 
gested that the hole conductivity is given by a thermally 


In the present work we have studied the oxidation and 
hydration of an acceptor-doped proton-conducting per¬ 
ovskite oxide, BaZrOg, in contact with water vapor and 
oxygen gas. Charge carrier concentrations have been de¬ 
termined for different temperatures and partial pressures 
based on data from first-principles modeling. 

Two different methods have been employed that im¬ 
prove upon the conventional PBE functional with regard 
to the description of band gap and band edges, namely 
the PBEO hybrid functional and PBE-l-GoWo calcula¬ 
tions rooted in many-body perturbation theory. 

We find that the hydration reaction is exothermic 
and well described by both PBE and PBEO. In¬ 
cluding the band edge shifts from GqWo calculations 
(PBE-|-x[GoWo]) does not change the energetics for the 
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hydration reaction. 

For the oxidation reaction, however, the different ap¬ 
proximations predict qualitatively different results. With 
PBE the reaction becomes exothermic while it is en¬ 
dothermic with PBEO and PBE-|-x[GoWo]. The exother¬ 
mic PBE behavior yields large hole concentrations when 
lowering the temperature even under hydrated conditions 
and the oxide can not become completely hydrated, in 
disagreement with experiments. For the exothermic na¬ 
ture of PBE to be consistent with the experimental data 
for the hole conductivity the hole mobility has to increase 
more rapidly than the decrease in hole concentration. 
Such a temperature dependent hole mobility is unlikely. 
We conclude that only the endothermic behavior with 
PBEO and PBE-l-GoWo can be made consistent with ex¬ 
perimental data of charge carrier concentrations and hole 
conductivities. 

In summary, PBE gives a good description for the hy¬ 
dration reaction but to model the oxidation reaction im¬ 
proved approximations have to be used. Here we show 
that the PBE-l-GoHo method and hybrid functionals are 
two viable alternatives and we present a theoretical ap¬ 
proach, which in a consistent way describes both hydra¬ 
tion and oxidation of proton conducting acceptor-doped 
perovskites. 
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Appendix: Band structure alignment 

To determine the shift of the VBM and CBM between 
PBE and PBEO the band structures need to be aligned. 
Such an alignment can be done with respect to a common 
reference potential, such as the vacuum level Wac'^ A 
schematic representation of the alignment is depicted in 
Fig. [TT] If the band structures are aligned with respect 
to this reference then the shifts of the VBM and CBM 
are given by the difference in the ionization potential IP 
and electron affinity EA, respectively, 

AevBM = - IpP®™ (A.l) 

AecBM = EA^®® - EA^®™, (A.2) 

where IP = Kac - cvbm and EA = Kac - ccbm- 

To determine Wac and consequently IP and EA a sur¬ 
face calculation has to be performed. Such a calculation 
requires a supercell containing a sufficiently long slab of 
BaZrOa so that the core of the slab becomes bulk-like, as 


PBE PBEO 



Figure 11. Schematic representation of the band structure 
alignment between PBE and PBEO. 


well as enough of vacuum, in order to reach the vacuum 
level. An important aspect of this approach is that the 
vacuum level of the slab system, Vvac,siab, is not the same 
as the desired vacuum level due to ionic and electronic 
relaxation at the surface of the slab and thus can not 
directly be used as vacuum level in the alignment pro¬ 
cedure. To obtain the actual vacuum level these surface 
contributions need to be removed: 

Kac = Kac.slab — AI4l — A Von, (A.3) 

where At4i and AVon are contributions from electronic 
and ionic relaxation at the surface, respectively. In the 
following only electronic relaxation is considered, hence 
A Von = 0. The desired IP can thus be extracted from 
the slab system according to 

IP = VvaCjSlab CyBM.slab AVl — IPsurf AVl (A. 4) 

and together with Eq. (jA.l|) we obtain the shift of the 
VBM according to 

Acvbm = IP™f-IPrnT-(AFer'=-AV™). (A.5) 

In the same manner we obtain the following expression 
for the CBM shift 

AecBM = EA^f - EA™ - (AV^"^ - AV,™™). 

The electronic relaxation at the surface gives rise to a 
surface dipole (see Fig. [ni. If we denote the difference 
in the planar averaged (in the xy-plane) charge density 
between the surface and bulk systems Ap{z), where 2 ; is 
the axis perpendicular to the surface, then the potential 
arising from the surface dipole can be calculated from the 
expression^i 

AVei = (A.7) 
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Figure 12. Schematic representation of how the surface dipole 
charge density Ap is obtained. The Ap-curve (green) has been 
multiplied with a factor of 5 for clarification. 


where eg is the vacuum permittivity, A is the unit area, 
and p is the electric dipole moment 

P= j Ap{z){z - zo)dz, (A.8) 

with Zq denoting the center of mass. 

In an actual calculation A.p(z) is obtained in the follow¬ 
ing manner (for a schematic representation see Fig. 1121) . 
First, a 1 X 1 X n supercell is constructed and the cor¬ 
responding charge density Pbuik(-2) is determined. Half 
of the atoms are then removed resulting in a supercell 
containing a 1 x 1 x ^ slab and an equal amount of vac¬ 
uum. Subsequent electronic relaxation yields the charge 
density psurf(-z)- Ai.p(z) is then obtained as the difference 
between these charge densities. 


where Pbuik('Z) has been truncated and set to zero at the 
same position as the surface in the slab supercell. 

To determine the VBM and CBM shifts we have con¬ 
sidered the [001] surface with both Zr02 and BaO ter¬ 
minations. We have used n = 9, which corresponds to 
a slab consisting of four and a half unit cells, where 
both surfaces (the second surface arise from the peri¬ 
odic boundary conditions) have the same termination. 
We use the same computational setup as described in 
Section imi however, only one /c-point is used in the z- 
direction. A summary of the results is given in Table IVl 
Using Eq. (IA.5I) we obtain VBM shifts of —1.44eV and 
— 1.38eV for the Zr02 and BaO-terminated surfaces, re¬ 
spectively. These shifts are in very good agreement with 
the VBM shift of —1.42 eV obtained by directly compar¬ 
ing PBE and PBEO results for the bulk. For the CBM 
shifts we obtain 0.83 eV for the Zr02-terminated surface, 
which compares well with the direct value of 0.80 eV. For 
the BaO-terminated surface, however, the CBM shift is 
only 0.34 eV. This discrepancy is likely related to the 
fact that the conduction band consists of zirconium d- 
states (see Fig. ED, which are not present in the surface 
layer for the BaO termination. In all, the results ob¬ 
tained here demonstrate the proper alignment of PBE 
and PBEO band structures (at identical lattice constant 
and using the same pseudopotentials) shown in Fig. [2] 

Table V. Difference between PBE and PBEO results for the 
bulk system as well as both terminations of the [001] surface. 
Equivalent band edge shifts for the different systems are given 
in bold. Energies are given in units of eV. 


Quantity 

Bulk 

Zr02 

BaO 

AevBM 

- 1.42 

-1.37 

-1.26 

AecBM 

0.80 

0.90 

0.36 

AEgap 

2.22 

2.26 

1.62 

tt-PBE tt-PBEO 

^ vac ^ vac 


0.10 

0.14 

-(akt®"= - AKr®°) 


-0.17 

-0.27 

AevBM from Eq. (IA.5II 


-1.44 

-1.38 

AecBM from Eq. (lA.6|l 


0.83 

0.34 


— Psurf(-2^) Pbulk(-2^); (-^■9) 
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Supplemental Material: Implications of the band gap problem on oxidation and 
hydration in acceptor-doped barium zirconate 

Anders Lindman, Paul Erhart, and Goran Wahnstrdm 
Department of Applied Physics, Chalmers University of Technology, SE-flS 96 Gothenburg, Sweden 

In this supplemental material we derive explicit expressions for the equilibrium constants of the 
hydration and oxidation reactions in acceptor-doped perovskite oxide materials. These expressions 
are used in the main article to make connections between the theoretical modeling and experiments. 

For the oxidation reaction we consider two limiting cases: high and low hole concentrations. 


I. OXIDATION REACTION 

Consider first the oxidation reaction 

1 02 (g)+ v^*^2h- + Og. (1) 

This reaction is associated with the Gibbs free energy of 
formation 


AGox(T,poJ = 2Me(T)-AGv, (2) 

where iie{T) is the electron chemical potential with re¬ 
spect to the valance band edge and AGv is the Gibbs 
free energy of formation for a doubly charged oxygen va¬ 
cancy. The corresponding equilibrium constant is given 

by 


ifox(T) 


\ - 1/2 2 
PO 2 \ 

PO 2 ) Cv 


(3) 


where n\,, c„ and cq denote hole, oxygen vacancy and 
oxygen ion concentrations, respectively, and po^ (Po^) i® 
the (reference) oxygen partial pressure. 

In the dilute limit the concentration of an arbitrary 
defect is given by 

(4) 

where Adef is the number of defect sites in the unit cell 
with volume I 4 and AGdef is the Gibbs free energy of 
defect formation. For the oxygen vacancy concentration 
we get 


Cy 


^ AG°AT)/kT -2ne{T)/kT (PoA 


(5) 


where 


AGAT) = AHAT) - TASAT) ( 6 ) 


is the standard Gibbs free energy of formation for the 
oxidation reaction. The site restriction on the oxygen 
sublattice implies that 


= 


and consistent with the dilute limit we have cq = N^/Vc, 
hence 

KAT) = n2e2Me(T)/fcTg-AG“JT)/fcT^ (g) 

The electron chemical potential is determined from the 
charge neutrality condition, which in general can be writ¬ 
ten as 

^ ^ ^defCdef(/le) ll'e(/le) ll'h(/le) — 0; (9) 

def 

where Ue is the electron concentration. In this case we 
have 


2cv “t" ri\, — ca, (^b) 

where ca = xxjVc is the acceptor dopant concentration 
and xa is the fractional acceptor dopant concentration. 

To obtain an explicit expressions for the equilibrium 
constant a model for the density of states (DOS) g{e) 
has to be introduced. We assume a simple parabolic ex¬ 
pression for the valence band 


5(e) 


5o-\/evBM — e, e < evBM 
0, e > evBM 


( 11 ) 


with 




( 12 ) 


where is an effective mass for the hole. 

We will now consider two limiting situations; low and 
high hole concentrations. 


A. Low hole concentration 

We first consider the situation where the hole concen¬ 
tration is low: nil “C u-v This is the situation when the 
PBE-|-x[GoWo] and PBEO approximations are used. In 
this case the electron chemical potential will be located in 
the band gap, pe > 0, and we will assume that pe kT. 
In that case the hole concentration can be obtained using 
the Boltzmann approximation 

/ evBM _ 

5o\/evBM - e[l - /(e,Me)] de 

-CO 


CQ + C, 


( 7 ) 


( 13 ) 
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where nvB(T) is the effective state density for the valence 
band, which is given by 


nvB{T) = 2 


V 2Trh^ ) 


(14) 


With nil “C Cv the charge neutrality condition implies 
that 2cv = CA- With Eq. dH) we then get the following 
expression for the electron chemical potential 


„2Me(T)/feT _ 


PO 2 

Po, 


- 1/2 


2W 


XA 


:!_^AG°JT)/kT 


(15) 


and the hole concentration is thereby given by 


nh{T,po2) 




nvB{T)e 


-AG°^{T)/2kT 


(16) 


Finally, for the equilibrium constant we obtain 


With nil ^ Cv the charge neutrality condition implies 
n-h = ca and using Eq. (1^ we get after some algebra the 
following expression for the electron chemical potential 


i^e{T) = a 



(24) 


with a = (3cA/2go)^'^^- This further implies that 

K,^{T) = cle-^^(^oAT)+2MT)]/kT ( 25 ) 

and 


= AH^^iT) + 2 a 




■ (26) 


The hole concentration gives rise to an effective oxida¬ 
tion enthalpy that is larger than the standard oxidation 
enthalpy AH°^{T). 


K^^iT) = [ny^{T)fe-^^°^AT)/kT^ASl^(T),k^ (^ 7 ) 

If we define an effective oxidation enthalpy 

AH:^^\T) = (18) 

we get from Eq. (II3 that 

Ai/r®(T) = AH:^{T) + 3fcT, (19) 

where we have made use of the thermodynamic relation 

^AH:^{T) = T-^AS:^{T). ( 20 ) 

In the same way, if we write 

nil = (21) 


II. HYDRATION REACTION 

Next we consider the hydration reaction 

H20(g)+v^* + 0g v^20H^. (27) 

The Gibbs free energy of formation for this reaction is 
given by 


AGiiydr(T,pH2o) = 2 AGh - AGv (28) 


where AGh is the Gibbs free energy of formation for 
a proton and PH 2 O is the water partial pressure. The 
equilibrium constant of Eq. is given by 


^hydr(T) 


f PH2o \ ^ Cb 

\PB 20 ) CvCo’ 


(29) 


for the hole concentration, we get an effective energy 


AEf{T) = -k 


din nil(T) 

d{llT) 


AH°^{T) 3fcT 
2 


( 22 ) 


B. High hole concentration 

Now we consider the case with high hole concentration: 
nil Cv This is the situation using the PBE approx¬ 
imation at low temperatures. In this case the electron 
chemical potential will be located in the valence band, 
/ie < 0. We define Ve = —pe and assume that Ve ^ kT. 
In this case the hole concentration can be obtained using 
the Sommerfeld expansion 


where ch is the proton concentration. 

The site restriction now also includes the proton con¬ 
centration 

Co + Cv + Ch = , (30) 

Cc 

however, in the dilute limit we obtain the same expression 
as for the oxidation reaction, cq = N^/Vc- Using Eq. (g]) 
for Cv and ch we obtain 

Aii,dr(T) = 

\PH2 0/ -'''v 

= (31) 


nh(l'e) 


5oV evBM 




3/2 


1-f 


- e [1 - /(e, Ve)] de 



(23) 


where 

AG;(,dr(r) = AH°^^^{T) - TAS°^^^,iT) (32) 

is the standard Gibbs free energy of formation for the 
hydration reaction. 



















